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ABSTRACT 



Complete solutions to the equations of motion are possible using 
numerical schemes. The model developed by Galt (1972) was investigated 
for stability and convergence. A series of runs indicated that the 
scheme, the three level Adams-Bashforth method for the new value of 
vorticity and then the successive over-relaxation method for the stream 
function, was stable for transport not exceeding 100 Sverdrups and a 
time step of one-half pendulum day. Convergence was determined by 
comparisons, with known analytic solutions, of phase velocity of Rossby 
waves, and for the steady state conditions, patterns of the stream 
function and the vorticity. Improved accuracy of the model was attained 
by using a finer grid. 
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I. INTRODUCTION 



The study of the dynamics of oceanic flow has been attempted by 
using the equations of motion. Since solutions to those complex equa- 
tions in their complete form was not probable it became necessary to 
simplify them by making highly restrictive assumptions that allowed 
for some of the functional relationships to be defined as constants and 
the non-linear terms to be disregarded. Ekman (1905) developed the now 
classic boundary layer solution, known as the Ekman spiral, which 
showed the effects of wind stress acting on a free surface. Stommel 
(1948) was able to reduce the complex equations to a balance between 
wind stress, a parameterization of friction, and Coriolis force. His 
solution dramatized the phenomenon of western intensification. 

The derivation of the vorticity equation from the vertically inte- 
grated equations of motion .removed the necessity of defining the press- 
ure field at the surface and simplified the boundary conditions. Munk 
(1949) developed a solution to a simplified form of the vorticity equa- 
tion which contained no non-linear terms, was invariant with time, 
assumed a zonal wind stress, and assumed a constant eddy coefficient. 

In both solutions the features Stommel observed in his solution were 
evidenced and, additionally, Munk (1949) observed the counter current 
region as well . 

Development of analytic solutions to the equations of motion 
necessarily demanded regular, closed boundaries, and for some forms, a 
constant depth. In each solution mentioned, an ocean of constant depth 
was assumed or, as in Ekman's solution, the ocean was considered as 



7 



unbounded in depth. To study, analytically, ocean basins of irregular 
boundaries, with inflow and outflow, and variable depth lead to an 
impossible task. 

Early work in developing these solutions was usually attempted for 
ocean basins that were rectangular in shape. A need to investigate the 
effects of a laterally sloping boundary was recognized, as evidenced 
by the more recent work, by Munk and Carrier (1950) which was done 

on a triangular shaped basin. Studies of the effects of bottom topo- 
graphy on ocean circulation have been presented by Warren (1963) and, 
more recently, by Veronis (1966). The solutions Veronis developed 
were to the time dependent, non-advective equations with no friction or 
wind stress and fora rectangular basin. 

It was clearly indicated that with each solution a different pheno- 
menon of oceanic flow could be described using the vorticity equation 
or the equations of motion. The optimum stiuation, then, would be to 
obtain a complete solution to these equations so that instead of being 
restricted to solutions which describe only one of the phenomenon a 
complete interaction of advection, Coriolis force, friction, and wind 
stress could be studied as a function of time. 

Since a numerical solution to the equations of motion removed the 
necessity of simplification and allowed for a complete solution, various 
models have been developed. Usually the only simplifications involved 
were to specify the density field and friction as constant. More 
recently Galt (1972) described development of a numerical model of a 
constant density form of the vorticity equation. 
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The non-dimensional equations that were used are represented by 



equations (1) and (2). 
Vortici ty : 
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= 


vorticity 


h 


= 


depth 


f 


= 


Coriolis parameter 


K 


= 


lateral kinematic eddy coefficient 


T 


= 


wind stress 




= 


horizontal gradient operator stream function 



In this form the non-dimensional ized constants a and 6 were defined: 

a = '1 'q/FL^D 

6 = K/FL^ 

Where: 

= a representative value of the stream function 

F = a representative value of the Coriolis parameter 

D = mean depth 

L = the grid length 

The non-dimensional wind stress x was related to x', the dimensional 
form, by equation (3). 
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If 




Where: 



p = density 

The non-dimensional stream function was then defined: 

’I' = 'J'' (FL^D) = '1'^ 

Using the non-dimensional forms of the equations the model was 
constructed. The new value of vorticity was obtained by the three level 
Adams-Bashforth differencing method. 

The stream function was then relaxed to within acceptable limits by 
the method of successive over-relation. This scheme was applied to a 
polygonal computational molecule. 

This computational grid allowed for the specification of irregular 

boundaries and for inflow and outflow. Newtonian, lateral friction was 

2 

parameterized in terms of the horizontal Ekman number, K/FL . 

For models which have been developed it has been necessary to, first, 
determine their stability and second, to determine their convergence to 
the true solution. Since complete, analytic solutions were not available, 
these models were run for conditions which were representative of 
simplified equations to which solutions existed. Comparisons were then 
made to assess the validity of the models. 

For some finite difference schemes it is possible to mathematically 
solve for errors generated and determine rates of convergence as well as 
showing what bounds there must be on the spatial and temporal steps to 
ensure stability. In referring to the model Galt developed it was not 
possible to make general calculations for convergence and stability. 

As indicated by Baer and Simons (1970) a method to investigate numerical 
stability in a model can be effected by monitoring a conservative 
property such as energy and total vorticity. The method used to determine 
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convergence was by a comparison of the patterns of vorticity and the 
stream function with known analytical solutions. This method, for 
example, was used by Blandford (1970). 

Establishment of the characteristics of Galt's model in terms of 
stability and convergence was done in two parts; 1): For stability, a 

form of kinetic energy, total vorticity, and the maximum value of the 
stream function were monitored; 2): To check convergence a comparison 

was made between known analytic solutions and the solutions predicted 
by the model. These comparisons used patterns of the vorticity and the 
stream function, and of the phase velocity of Rossby waves. These 
comparisons were made with two different model configurations. 

Once the basic shape and boundary conditions were determined the 

initial input was either constant or variable depth, distribution of 

the Coriolis parameter, magnitude and direction of the wind stress, 

stream function and vorticity patterns, and the local rate of change and 

vorticity one time step back. The parametric specification of density 

and frictional coefficients, which were assumed constant, and the 

physical dimensions of a particular configuration were determined by 

7 2 

the magnitudes of the non-dimensional izing constants , 'I'^/FL D and K/FL . 
Additionally, for those runs when these constants were zero the 
dimensions of the model were determined by the distribution of the 
Coliolis parameter, magnitude of the wind stress, and bottom topography. 

The first configuration was a rectangle shape that could only be 
approximated by the model due to the use of the polygonal molecules 
for the computational grid. The boundary conditions were set up so as 
to represent an enclosed basin with free slip boundaries with no 
friction in the model. The non-linear terms were not used and the wind 
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stress set to zero so that this solution could be compared to the free 
wave solutions of Veronis (1966). One run was with constant depth and 
a variable Coriolis parameter and the other run was with variable 
depth and a constant Coriolis parameter. 

The second series of runs were done using a parallelogram shape 
which had closed regular boundaries. The east-west boundaries were 
inclined sixty degrees from the horizontal, measured in a positive 
sense. This was run with constant depth and a variable Coriolis 
parameter. The north-south boundaries were set to free slip conditions 
and the east-west boundaries were set to no slip conditions with 
friction in the interior. A zonal wind stress set in as a consine func- 
tion with maximum values at the north-south boundaries. The model did 
not accurately convert this to the curl of wind stress and, therefore, 
all successive runs were made with a direct input of the curl of the 
wind stress. These particular requirements were so that a comparison 
with a modified solution of Munk and Carrier's (1950) could be made. 

The investigation of the effect of reducing the spatial step size 
was done for both configurations by increasing the number of points 
and appropriately scaling the input parameters. Reduction of the grid 
length to one-half was accomplished so as to better define the 
velocities and western boundary layer. 
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II. STABILITY ANALYSIS 



Baer and Simons (1970) indicated that for the Adams-Bashforth 
iteration scheme, when used on a non-linear equation, it was sufficient 
to determine stability for a linear form of that equation since the on- 
set of instability appears to be the breakdown of the linear solutions. 
Since determination of computational stability was desired it was suffi- 
cient, (Baer and Simons [1970]), to monitor the conservative properties 
which were, kinetic energy, total vorticity, and the maximum value of 
the stream function. For the stable case. Figure 1, these properties 
were bounded over the period of integration. Figure 2, the unstable 
case, indicates some of these values as growing unbounded after a certain 
integration time. The model's response for the stable situation was a 
comparatively short spin-up time, as indicated by the values of the 
properties monitored. The -kinetic energy and maximum value of the 
stream function reached a maximum value at this time and were then 
observed to oscillate throughout the remainder of the integration time. 
This response was similar to the results reported by Blandford (1970) 
for a model in which was specified a lateral viscosity and no slip 
boundary conditions. Total vorticity reached as maximum value later in 
the integration then remained constant throughout the remainder of the 
integration time. For the unstable situation the kinetic energy after 
a period of time began to grow in an exponential manner with oscilla- 
tions of the stream function increasing. The vorticity, for this 
situation, grew steadily throughout the integration time. 
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NON-DIMENSIONAL MODEL TIME 



To more accurately determine the limits of stability various 
conditions were imposed on a constant depth model basin with basic 
configuration as shown in Figure 3. 





Figure 3. Basic Model Configuration for Stability Analysis 

For all but two runs mixed boundary conditions were used with the north 
and south boundaries being free slip. Of the other two runs, one was 
run with free slip boundary conditions, the other with no slip boundary 
conditions . 

The distribution of the wind stress over the model basin was zonal 
in nature. The direction of the wind stress was the cosine function in 
the north-south direction with the maximum stress in the westward 
direction at the south boundary and the maximum eastward stress at the 
north boundary. 

The first series of runs were made with the constant ('I’q/FL D) set 
to zero v/hich resulted in the linear form of the equations. The wind 
stress was set at a representative value and the Coriolis parameter was 
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allowed to vary over the constant depth model basin. The physical dim- 
ensions of the basin were determined from this distribution. Table I 
shows the stability for a range of friction, K. For the grid length 

L = 4.6 X lO^m, F = lO’^sec"^; K ranges from approximately 7.5 x 10^ to 
7 2-1 

7.5 X 10 cm sec . The effect of varying the magnitude of the wind stress 
for this test had the effect of varying the transport, and by using 
large values for. the wind stress the model was unstable. 

Table I. Stability Dependence on Friction and Time Step. 







DT 


KDT/FL'^ 


* 


STABLE 


2.86 X 


10-3 


0.5 


1 .43 


X 


10"3 


Yes 


2.86 X 


10"^ 


0.5 


1.43 


X 


10-3 


Yes 


2.86 X 


10"^ 


1.0 


2.86 


X 


10-^ 


Yes 


2.86 X 


10"^ 


2.0 


5.72 


X 


10-^ 


Yes 


2.86 X 


10"^ 


3.0 


8.58 


X 


10-^ 


? 


2.86 X 


10"^ 


5.0 


1 .43 


X 


10-3 


? 


2.86 X 


10"^ 


7.5 


1.71 


X 


10-3 


No 


2.86 X 


10"^ 


10.0 


2.86 


X 


10"3 


No 



Op 1 

Allowing K to remain constant at 7.5 x 10 cm sec" the wind was 
varied as shown in Table II. 

Table II. Stability Dependence on Transport 



curl DT STABLE 



-sin(y) 


0.5 


No 


0.0 


0.5 


Yes 


-lO-^sin(y) 


0.5 


Yes 


-lO-^sin(y) 


0.5 


Yes 


-10-''sin{y) 


. 0.5 


Yes 
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By using the non-dimensional izing coefficient for the wind stress, 




where p is density, F a representative value of the Coriolis parameter, 

O 

L the grid length, in Sverdrups (10 m /sec) was determined. The 
maximum absolute value of the non-dimensional wind stress, t, was 1 .0 . 
with p = 10"^ grams m"^, F = 10"^sec"\ L = 4.6 x lO^m the runs depicted 

5 

in Table II represent a range of maximum transport from zero to 10 
Sverdrups . 

The effects of non-linearities on the stability of the model are 

7 2 

shown in Table III. Keeping K constant at 7.5 x 10 , 4 'q/FL D was allowed 
to vary with varying wind stress. 



Table III. 


Stability Dependence 
Time Step. 


on Non-linear 


Terms and 


curl T 

A 


'I'q/FL^D 


DT 


STABLE 


-10"^sin(y) 


10'^ 


1.0 


Yes 


0.0 


10'^ 


1.0 


Yes 


-10"^sin(y) 


10'^ 


1.0 


No 


-10"^sin(y) 


10'^ 


1.0 


Yes 


-lO'^in(y) 


10'^ 


5.0 




-10‘^sin(y) 


10'^ 


1.0 


Yes 


-10"^sin(y) 


10-3 


10.0 


No 



2 

Inclusion of the non-linear terms by ’I'q/FL D where D is the mean model 

-2 2 

depth, ranged from 10 to 10 Sverdrups. This also scaled the 

• -13 

magnitude of the wind stress and it ranged from approximately 2 x 10 

-9 2 -2 

to 2 X 10 m sec for its maximum values. Similar, unstable, response 

2 

for large transport was obtained by setting 'J'q/FL D to 10.0 which 

5 

represented approximately 10 Sverdrups. 
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Table IV indicates the final three runs made for determination of 
computational stability. Again, using the linear form of the equations 

o O 

and K approximately 10 cm sec , the maximum transport for the wind 
stress used was approximately 10 Sverdrups. 



Table IV. Stability Dependence on Time Step. 



Curl 


K/FL^ 


DT 


STABLE 


-10”^sin(y) 


5 X 10"^ 


1.0 


Yes 


-10"^sin(y) 


5 X 10"^ 


5.0 


? 


-10”^sin(y) 


5 X 10'^ 


10.0 


No 



An attempt to compute computational stability was made by using; 

fc ^ + K -^)< 1.0 

\ (AX)^/ 

The term, C — where C is a characteristic velocity, was not important 

AX 

in the series of linear runs. The term K was computed as shown 

in Table I with the results of a somewhat larger bounds for At. For 

the transport range of 10 to 100 Sverdrups and a coefficient of friction 

8 2-1 4 

on the order of 10 cm sec a time step of 10 seconds ensured stability. 

This represented a one-half pendulum day. The effects of the different 
boundary conditions mentioned previously was to change only the model 
spin up time. This observation was similar to that reported by Blandford 
(1970). 
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III. ROSSBY WAVE PROPAGATION 



Convergence of the model to a steady state analytic solution was 
determined by allowing the model to generate time dependent, free wave 
solutions and allowing the transients to die out leaving only the steady 
solution. These Rossby waves can occur in an ocean basin and are of 
two types. The classical Rossby wave can exist in a constant depth 
ocean due to the change in the Coriolis parameter with latitude. Con- 
sidering the situation where there is no friction and no wind stress, 
with the Coriolis parameter increasing from south to north, a Rossby 
wave will tend to propagate along lines of a constant Coriolis parameter. 
For this case, then, the wave would propagate in an east to west direction. 

The second type of Rossby wave that can exist in an ocean basin is 
the topographic Rossby wave. Consider, now, a constant Coriolis para- 
meter but with a simply varying depth which increases from east to 
west. The topographic Rossby wave will propagate as before except now 
along lines of constant depth in an east-west direction. For a simply 
varying depth, increasing from east to west the topographic Rossby 
waves will propagate from south to north across the ocean basin. These 
waves can be described by solutions to the potential vorticity equation. 



Considering that potential vorticity is conserved following a parti- 
cle of water’, equation (4) states that for free wave solutions the second 
term is zero, tiat is, f/h is constant. Using equation (4), Veronis 
(1966) developed a time dependent, harmonic, constant depth solution 




C = vorticity 



(4) 
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for the stream function, ip (x,y,t). 



<P = 'pQ sin( Ax + cut) sinax sinyy (Veronis[1966]) (5) 



Equation (5) describes sinusoidal waves propagating across a basin in 
the -X direction enclosed in a sine envelope. For the first mode in the 
X and y directions the following constants were defined: 

A = B/2 

a = tt/L 

Y = ir/M 1 

» = e/2[J + 

af 

e = constant = — 

oy 

L = basin dimension in x direction 

M = basin dimension in y direction 



For use in initializing the model non-dimensional ization of equation (5) 
was necessary. By non-dimensional izing L, M and a with a grid length, 
aL, and to with F, the following non-dimensional constants were defined: 



A 

I 

a 

I 

Y 



r '2 + 

La +Y 



2 



1 



(tt/L) aL 



(tt/0.866L) aL 



I I 

to = BaL/2a F 



aL = grid length 

L = basin dimension 

F = a representative value for the Coriolis parameter. 



The phase velocity for this solution is to/A which means that the non- 
dimensional phase velocity is to'/A' = 3 aL/2x^ F, or BAL/2F[a ^ + y ^]- 
» A .second solution developed by Veronis (1966) was for a linearized 
form of equation (4), but with a constant Coriolis parameter and with 
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• kx 

variable depth that was assumed to be h = h^e . This solution, 
equation (6), described sinusoidal wave patterns: 

= i>o sin(xy + wt) sinax sinyy (6) 

propagating in the positive y direction which are enclosed in a sine 
envelope. The envelope is increasing in an exponential manner in the 
X direction. Again, for the first mode, the dimensional parameters were 



defined. 










X 


= 


fk/2o3 




a. 


= 


ir/L 




Y 


= 


ir/M ^ 




0) 


= 


fk/(4[a^ + Y^] + 




f 


= 


constant Coriolis parameter 




L 


= 


basin dimension in x direction 




M 


= 


basin dimension in y direction 



Using aL to non-dimensionalize L, M, and A, and to with F the following 
parameters were defined: 

k' = kAL ^ 

co' = k'/{4[a ^ +Y ^] + k'^)^ 

a = (ir/L) AL 

Y = (ir/0.866L) aL 

I I 

The non-dimensional phase velocity, w /A , was obtained: 

co'/a' = 2k'/{4[a'^ +y'^] + k'^)^ 

The rectangle shaped basin which was desired could not be exactly 
specified as depicted in Figure 4. This model basin shape had sawtoothed 
east-west boundaries due to the computational grid being made up of 
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polygonal molecules. Scaling the model basin to physical dimensions was 
done in the following manner. For the constant depth configuration the 
model was scaled on the distribution of the Coriolis parameter. 




AM 



8 = — = n = 9 19 

^ ay nAM » 

aF = 0.75 X 10"^sec“' 
a = 60° 



using a representative value for $ = 
1.8 X 10‘^^ 'lsec"\ 



am = 



aN = 



.75 



(9)(1.8xl0'^b 



= 4.63 X 10 meters 



aM 



sin 60' 



= 5.35 X 10 meters 



The configuration for which the Coriolis parameter was constant and 
the depth varied as in equation (6) determined the physical dimensions 
of the model in the following manner; 

L _ l_ 

h = h„e 
0 



h 

In 



X = 



= nAN , n = 9,19 

k = 2.2 X 10“V^ 
h^ = 3.0 X lO^m 
h = 1.0 X lO-^m 



an = 5.56 X lO^'m 

am = sin 60 °(aN) = 4.83 x lO^m 



Since free wave solutions were desired, all of the boundary condi- 

. 2 ? 
tions were set at free slip. The constants K/FL and 'f’^/FL D were set 

to zero so that the model had no friction and so the linear equations 

could be solved. Since the wind stress was also zero for these runs 

initialization of the stream function and corticity patterns was 



23 



necessary. To accomplish this, the non-dimensional i zed solutions 
developed earlier from equations (5) and (6) were solved at the non- 
dimensional time 5 and the resulting patterns used as input to initialize 
the model . 

In general , the model propagated the two types of Rossby waves across 
the basin with the kinetic energy and total vorticity remaining essen- 
tially constant. The fluctuations in kinetic energy that were observed 
had their maximum value when the wave, depicted in Figure 4, reached 
the center of the basin. This was expected since the waves were pro- 
pagating in a sine envelope which reached its maximum at the center of 
the basin. The influence of the boundaries on the waves appeared to 
cause them to disperse laterally. This observed influence quickly 
disappeared as the waves moved to the interior of the model 




< Direction of Movement 



Figure 4. Example of Rossby Wave Propagation. 
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Reducing the grid length to one-half its former value by doubling the 
number of points on a side did not change these observed characteristics. 
Subsequent analysis indicated that these waves were more accurately 
defined with the decreased grid length. 

The phase velocity for both series of runs was determined by using 
Figures 5 through 8. The first comparisons with analytic solutions are 
shown in Table V. This was for the configuration denoting a basin of 
constant depth and with variable Coriolis parameter. 



Table V. Comparison of Rossby Wave Velocity. 



NUMBER OF GRID GRID LENGTHS/10 NON-DIMENSIONAL PHASE Vel 
POINTS/SIDE TIME STEPS Model Analytic Solution 



10 



1.2 



.12 .169 



20 



2.4 



.24 .356 



The percent error for ten grid points was 29 % and that for twenty grid 
points was 17 %. 

The configuration of a basin with an exponentially varying depth 
and with a constant Coriolis parameter was analyzed in a similar manner 
for topographic Rossby waves. Table VI indicates the comparison made 
with the analytic solution. 



Table VI. Comparison of Topographic Rossby Wave Velocity. 



NUMBER OF GRID 


GRID LENGTHS/10 


NON-DIMENSIONAL 


PHASE Vel 


POINTS/SIDE 


TIME STEPS 


Model Analytic 


Solution 


10 


1.2 


.12 


.212 


20 


3.5 


.35 


.445 
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Figure 5. Non-dimensional Stream Function at Every Ten Time Steps for 
Phase Velocity Determination, the Classic Rossby Wave. 
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Figure 8. Non-dimensional Stream Function at Every Ten Time Steps for 
Phase Velocity Determination, the Topographic Rossby Wave. 
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The same trend was observed in that the percent error for ten grid 
points was 33 % and that for twenty grid points was 21 %. 

This series of runs established a very important trend for the 
model. Increasing the number of computational points increases the 
accuracy of the model for the free wave solutions. In approximate 
terms, doubling the number of computational points very nearly halves 
the percent error as compared to the known analytic solutions. This 
determination was subject to graphical interpretation which introduced 
possible errors. In an attempt to clarify the picture, the figures 
were normalized with respect to the sinax envelope but no significant 
improvement resulted. 



30 



IV. THE STEADY STATE SOLUTION 



The assumption that steady state conditions exist in an ocean basin 
was made by Munk and Carrier (1950) when they derived a beta plane solu- 
tion on a triangular basin of constant depth. This solution to the 
vorticity equation balanced frictional effects, Coriolis force, and an 
assumed wind stress distribution. To make this compatable with the 
configuration of the numerical model. Figure 9, the dependence of the 
basin width in the x direction on y was removed so as to make it a 
constant width. The solution was then corrected for the boundary 
conditions which were no flow normal to the boundaries. This solution. 



IJ 

h 




Wind Stress 



Figure 9. Basic Model Configuration 



in non-dimensional form is shown in equation (4); 
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with the following parameters defined: 
e = y/P 

p = (1-$^ tan^e)'^^^ 

0 = 60° 

$ = l-l/3((j)- <j)^) tan<}) 

<t>^ = latitude at which y = 0 

^0 

<|) = latitude at which y = tt 

a = x/y - (y/Y) tane 

Ya = width of the basin 

The modified solution is shown in Figures 10 and 11. These patterns 

were then compared to the steady state reached by the model . To 

parallel the model with this solution the boundary conditions were 

set at free slip conditions for the north and south boundaries and no 

2 

slip for the east and west boundaries. The constant, K/FL , curl of 
the wind stress, and the distribution of the Coriolis parameter were 
scaled with respect to each other so as to represent a constant depth 
basin with approximate physical dimensions of 5 x 10 km in the x 

3 

direction and 4.5 x 10 km in the y direction, with the curl of the wind 

-4 

stress being distributed as 10 sin(y) over the basin. The representa- 

8 2 

tive value for friction, K, using these values, was 10 cm /sec. 

As indicated by Figure 12 this configuration reached steady state 
at about 1 x 10^ seconds. This time is when the maximum value of the 
stream function was essentially at its steady state value with the vor- 
ticity and kinetic energy within 95 % of their steady state values. To 
determine the convergence of this to the analytic solution, patterns of the 
vorticity and of the stream function were drawn. Figures 13 and 14. The 
difference between these steady state patterns that the model predicted 
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Figure 10. Contoured Values of the Non-dimensional Stream Function 
for the Modified Solution. 



1 










34 



Figure 11. Contoured Values of the Non-dimensional Vorticity for 
the Modified Solution. 
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Figure 12. Model Response. 
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Figure 13. Steady State Pattern of the Non-dimensional Stream Function 
Predicted by the Model. 
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Figure 14. Steady State Pattern of the Non-Dimensional Vorticity 
Predicted by the Model . 
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and that of Munk and Carrier could probably be attributed to the differ- 
ence between the physical conditions which brought the two solutions to 
a steady state. A numerical model which Gates (1970) used to solve a 
slightly different form of the equations of motion produced similar steady 
state patterns. He attributed this to Rossby wave reflections from the 
boundaries with the most important reflection being the one from the 
waves which propagated- from east to west. These waves diminished in 
magnitude with each crossing as well as diminishing in magnitude as they 
propagated through the model. The basic differences, not considering 
the finite difference methods, between Gates' model and Galt's model 
were: all boundaries no slip as opposed to mixed free slip/no slip 
boundaries, and free surface as opposed to a rigid lid. These differ- 
ences did not cause the two models to vary to any great extent other 
than a longer spin up time and a phase shift in Galt's model. The most 
important parallel is that of the type of friction which was lateral 
and did not include any bottom friction. 

Blandford (1970), after considering different combinations of 
boundary conditions, was able to show by using either free slip or no 
slip lateral boundaries that bottom friction, expressed as an exponen- 
tial boundary layer, dissipated more energy than did lateral friction. 

By including this boundary layer in the solution Blandford noted that 
the Rossby waves which were generated with model spin up tended to be 
damped out so that they did not dominate the interior flow. The boundary 
layer also caused a reduction in the time at which the model reached 
steady state conditions. 

The response of this model, early in the integration, was to develop 
maximum transport near the center of the basin. This cell soon propagated 
to the western edge where it dominated the flow as the boundary layer. 



38 



After the boundary layer was established, Rossby waves began to initial- 
ize in the southeast corner and propagate to near the western boundary. 
These inertial waves remained near the southern boundary and became 
established in the flow as the weaker counter current region. Gates 
reported observations on these waves obtaining the phase velocity of 
5.30 meters/sec for their westward component. These waves, depicted in 
Figures 15 and 16, had as their western component of the phase velocity 
5.34 meters/sec for the larger grid and 4.83 meters/sec for the smaller 
grid. 

To better define the phenomenon observed in Galt's model, a reduc- 
tion of the spatial step size was attempted. Already determined, in 
Section II, was the improved accuracy of the solution with a smaller 
grid length. The response of the model for the stream function for 
different grid lengths, aL, was observed in Figure 17 as being a higher 
maximum value for (1/2 )aL. This difference was not observed in the 
comparison of kinetic energy. Figure 18. The difference, in the stream 
function, between lAL and (1/2)AL was attributed to the influence of 
transients for the larger grid length. A profile of the stream func- 
tion for the two grid lengths. Figure 19, indicated that for the larger 
grid length, lAL, the solution oscillates eastward of the counter 
current region. The smaller grid length, (1/2)AL, had a well defined 
boundary current and the weaker counter current region. The transport 
decreased rapidly eastward of this region with relatively little 
oscillation. Comparing this profile with the profile of the modified 
Munk and Carrier solution. Figure 20, it was observed that the model 
did in fact compare favorably with the analytic solution in its pre- 
diction of the western boundary current and the decrease in transport 
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Figure 15. Transient Rossby Wave. 
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Figure 16. Transient Rossby Wave, 
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*150 300 450 600 

NON-DIMENSIONAL MODEL TIME 

Figure 17. Comparison of the Non-dimensional Stream Function for IaL and 1/2 aL 
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eastv/ard of the counter current region. The contoured values of the 
stream function. Figure 21, supported this observation. 

The steady state solution generated by the model could be used to 
acequately indicate oceanic transport. The accuracy of this solution 
was dependent on grid length, the more accurate solution being obtained 
using the smaller grid length. 
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Figure 21. Dtstribution of the Non-dimensional Stream 
Function for 1/2 aL. 
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V. CONCLUSIONS 



The model was sufficiently stable for the lengths of integration 
time used. The model had to be scaled to reflect a transport from zero 
to 100 Sverdrups with a time step which represented a one-half pendulum 
day to ensure stability. To reach steady state conditions it was neces- 
sary to depend on this factor so that the transients would be damped out. 

The curl operation was not accurate enough to produce the correct 
results when operating on the wind stress. 

Rossby waves, either topographic or inertial, were satisfactorily 
described by the model. The phase velocity was more accurately defined 
with a refined grid. Although these results were dependent on graphical 
interpretation this observed trend persisted through the following 
series of runs in Section III. The errors were significant enough to 
warrant future investigations. 

The steady state solution generated by the model compared favorably 
with the modified Munk and Carrier solution. To reach this solution a 
long integration time was necessary. The spin up time was influenced 
by mixed boundary conditions which caused a longer time to reach steady 
state conditions. The transient solution the model developed was 
similar to that which was reported by Gates. This dominating influence 
could be more accurately described by including some form of bottom 
friction. These transients became imbedded in the flow and were signi- 
ficant in establishing the counter current region at steady state 
conditions. The oscillations observed in the stream function profile 
for larger grid length were due to these transients dominating the 
steady state flow. 
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APPENDIX THE FINITE DIFFERENCE EQUATIONS 



The specific details of the derivation of the finite difference 
equations that are given below can be found in Galt (1972). This deriva- 
tion incompasses a combination of work done by Winslow (1967), Sadourny, 
Arakawa, and Mintz (1967), and Williamson (1969). By considering a 
computational molecule. Figure 22, the solution for each term in equation 
(8) was transformed from a surface integration to a line integration by 
using Stokes theorem. 
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. -h vv„ 5 + ''h (f) 

( 2 ) 



( 8 ) 



(1) (2) (3) 

The value of term (1) at point L was considered as an average of the 
potential vorticity over the molecule times the value of the stream 
function; 

potential vorticity (PVO = 
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The value of term (2) at point L was considered to be the average 



value over the inner shaded area of Figure 22. 



67^5 = 6 



5L+I ^ jL-NH ^ jL-N ^ ^ 5L-I 



5L+N-I , 



Term (3) considered the average values of magnitude and direction of t 
along each side of the molecule, then stored as a constant array since 
it was invariant with time. 




The successive over-relaxation procedure was used for equation (2) 
where the residual was computed by equation (9) 




Res was determined by an average value of ^ '!') over the inner 

hexagon, Figure 22. 
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The amount of over-relaxation was determined by equation OO). 
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Figure 22. The Computational Molecule. 
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